#This to obtain:
#Figure_5a
#Figure_5b
#Figure_5c

library("ggplot2")
library("glue")
library('gridExtra')
library('grid')
library('glue')
library('MatchIt')
library('forcats')

setwd("/Users/bgpopescu/")
#setwd("C:/Users/bogdanp/")

#Reading Polygons
x<-st_layers(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb")


HRV_adm0 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm0_wgs")

HRV_adm3 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm3")
HRV_adm3<-st_simplify(HRV_adm3, dTolerance = 0.005)
HRV_adm3<-subset(HRV_adm3, select = c(ID_2))
HRV_adm3_data = read_excel("./Dropbox/Legacies_Project/Analysis/data/merge.xlsx", sheet=1, col_names = TRUE, skip = 0)

final_sp<-left_join(HRV_adm3, HRV_adm3_data, by = c("ID_2"="ID_2"))
data<-subset(final_sp, bosnia_border_mun==0)
data$dist1 <- as.numeric( with (data,ifelse(data$treat==1, 1, -1)))
data$dist2 <-data$dist1*data$krajna6_distance

data_selection<-subset(data, select=c(ID_2,
                               pct_military_1857,
                               railroads_1869_density,
                               pop_zadrugas,
                               zadrugas_exist,
                               dist2,
                               POINT_X,
                               POINT_Y,
                               bfe1, bfe2, bfe3, bfe4, bfe5, bfe6, bfe7,
                               treat,
                               pct_no_water,
                               share_lesseduc_25,
                               pov_rate_income))
  
trimmed <- data_selection %>%
  filter(!is.na(treat) &
             !is.na(POINT_X) &
             !is.na(POINT_X) &
             !is.na(pct_military_1857) &
           !is.na(railroads_1869_density) &
             !is.na(pop_zadrugas) &
           !is.na(zadrugas_exist) &
           !is.na(pct_no_water),
           !is.na(share_lesseduc_25),
           !is.na(pov_rate_income))

#Function to create graph for model
coeff_plot_models <- function(models, data_sample, y, model_labels, graph_title) {
  
  #Create empty variables
  confint1 <- rep(NA, length(models))
  confint2 <- rep(NA, length(models))
  estimate <- rep(NA, length(models))
  
  for (i in 1 : length(models)){
    specification = as.formula(paste(glue("treat~"), 
                                     paste(models[[i]], collapse = "+")))
    #Step 2: Matching procedure
    m.out1 <- matchit(specification, data_sample, method = "nearest", replace=F)
    
    #Step 3: Create a new DataFrame
    mmdata <- match.data(m.out1)
    
    t.val=t.test(mmdata[[y]][which(mmdata$treat==1)], 
                 mmdata[[y]][which(mmdata$treat==0)])
    confint1[i]=t.val$conf.int[1]
    confint2[i]=t.val$conf.int[2]
    estimate[i]=t.val$estimate[1]-t.val$estimate[2]
    
    t.val_unmatched=t.test(data_sample[[y]][which(data_sample$treat==1)], 
                           data_sample[[y]][which(data_sample$treat==0)])
    confint1_unmatched=t.val_unmatched$conf.int[1]
    confint2_unmatched=t.val_unmatched$conf.int[2]
    estimate_unmatched=t.val_unmatched$estimate[1]-t.val_unmatched$estimate[2]
  }
  
  df1 <- data.frame(model_names=model_labels, point_estimate = estimate, 
                   lower_b = confint1, upper_b = confint2)
  df2 <- data.frame(model_names="Cumulative \nEffect", point_estimate = estimate_unmatched, 
                    lower_b = confint1_unmatched, upper_b = confint2_unmatched)
  
  df <- rbind(df1, df2)
  df$order<-rownames(df)
  print(df)
  
  
  df_final<-df%>%
    mutate(name = fct_reorder(model_names, order))

  #return(paste(glue("df", paste(models[[i]], collapse = "_"))))
  actually_graph<- ggplot(df_final, aes(x= name, y=point_estimate))  + 
    geom_text(mapping = aes(label = round(point_estimate, 2)), nudge_x = 0.25,
              nudge_y = 0.25, size = 5, color ="red")+
    geom_point(size = 5)+
    geom_hline(yintercept = 0, colour = "black", size = 2, linetype = "dashed")+
    geom_errorbar(aes(ymin = lower_b, ymax = upper_b), size = 1, width = 0.1)+
    theme_bw() +
    labs(x = "Model", y = "Estimate")+
    #scale_x_continuous(name = "Model") +
    #scale_y_continuous(name = "Estimate")+ 
    theme(axis.text.x = element_text(size=14),
          axis.text.y = element_text(size=14),
          axis.title=element_text(size=14),
          plot.title = element_text(hjust = 0.5),
          legend.position = c(1, 0),
          #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
          #and c(1,1) corresponds to the "top right" position.
          legend.box.background = element_rect(fill='white'),
          legend.background = element_blank(),
          legend.text=element_text(size=12))+
    ggtitle(glue("{graph_title}"))
  print(actually_graph)
  ggsave(actually_graph, file = glue("Figure{figure_name}.jpg"), 
         height=10, width=20, 
         units = "cm", dpi = 300)
}

  
regressors<-list(c("pct_no_water", "Pct. Dwellings No Water, 2011", "_5a"),
                 c("share_lesseduc_25", "Pct. Less Educated People, 2011", "_5b"),
                 c("pov_rate_income", "Pct. People at Risk of Poverty, 2011", "_5c"))

specifications <- list(c('POINT_X', 'POINT_Y', 'zadrugas_exist', 'railroads_1869_density'),
                       c('POINT_X', 'POINT_Y', 'pct_military_1857', 'railroads_1869_density'),
                       c('POINT_X', 'POINT_Y', 'zadrugas_exist', 'pct_military_1857'))

names(specifications) <- c('eff_soldiers', 'eff_zadruga', 'eff_railroads')


specification_label<-c("Pct. Soldiers", 
                       "Zadrugas", 
                       "Historical \nUnderprovision \nGoods")

names(specification_label)<-c('eff_soldiers', 'eff_zadruga', 'eff_railroads')


setwd("./Dropbox/Legacies_Project/Paper/figures/")

for (y in regressors) {
  regressor_name  = y[[1]]
  regressor_label = y[[2]]
  figure_name = y[[3]]
  coeff_plot_models(models=specifications, data_sample=trimmed, y=regressor_name, 
                    model_labels=specification_label, graph_title=regressor_label)
}
